clc
clear
%%
%n>3
%DrawGrid(Grids,k) can draw cells in 2-d space.
%Continual learning classification method for time varying data
%Data:n-dimension and m types
%data:[d1,d2,...,dn,time,label],time is the data acquisition
%Grids:[d1,d2,...,dn,p,q1,t1,q2,t2,...,qm,tm]
%di:dimension
%qi:number of samples in one cell
%ti:the last active time
%data.time > grids.t learning. 
%All data.time is constant, CLCTVDAIS degenerate into a supervised learning classification method

%Load data. 
load DLD.mat

mingen = 4;  
dlttime = 20;
sym TrainSets;
sym TestSets;
sym CellSets;
sym CeSaSets;

%DLD
l1 = 100;
n = 1 : 2000;
n(1,1:l1) = l1;
D1 = [A';n;ones(1,n(end))]';                                                    %The lable of 1st class samples
D2 = [B';n;2.*ones(1,n(end))]';                                                 %The lable of 2ed class samples
D3 = [C';n;3.*ones(1,n(end))]';                                                 %The lable of 3rd class samples
tt1 = ceil(l1*0.8);                                                          %The number of training samples
tt2 = ceil(l1*0.8);
tt3 = ceil(l1*0.8);
di = 2;
t = 3;

for j = 1 : 3
    j    
    JG = [];
    t1 = randperm(l1);                                                     %Rrandomly selecte training data
    t2 = randperm(l1);
    t3 = randperm(l1);
    D11 = D1;
    D22 = D2;
    D33 = D3;
    Tr1 = D11(t1(1,1:tt1),:);                                               %The 1st class samples for training
    Tr2 = D22(t2(1,1:tt2),:);                                                %The 2ed class samples for training
    Tr3 = D33(t3(1,1:tt3),:);                                                %The 3rd class samples for training
    D11(t1(1,1:tt1),:) = [];
    D22(t2(1,1:tt2),:) = [];
    D33(t3(1,1:tt3),:) = [];
    Te1 = D11;                                                             %The 1st class samples for testing
    Te2 = D22;                                                             %The 2ed class samples for testing
    Te3 = D33;                                                             %The 3rd class samples for testing
    TrainData = [Tr1;Tr2;Tr3];
    TestData = [Te1;Te2;Te3];
    TestData = sortrows(TestData,di+1);
    
    if di > 3
        for ki = 1 : di-2
            TrainSets{ki} = TrainData(:,[ki,ki+1,ki+2,di+1]);
            TestSets{ki} = TestData(:,[ki,ki+1,ki+2,di+1]);
        end
        TrainSets{ki+1} = TrainData(:,[1,di-1,di,di+1]);
        TestSets{ki+1} = TestData(:,[1,di-1,di,di+1]);
        TrainSets{ki+2} = TrainData(:,[1,2,di,di+1]);
        TestSets{ki+2} = TestData(:,[1,2,di,di+1]);
        for ki = 1 : di
            [CellSets{ki} CeSaSets{ki}] = InitG(TrainSets{ki},mingen,t);
        end
        tt = size(TestSets{1},1);
        for ii = 1 : tt
            for ki = 1 : di
                %return affinity
                [CellSets{ki},CeSaSets{ki},JGI(ki,:)] = CTestfun(CellSets{ki},CeSaSets{ki},TestSets{ki}(ii,:),dlttime);
            end
            [Temp1,Temp2] = max(sum(JGI));
            for ki = 1 : di  %Adjust cell information
                [CellSets{ki},CeSaSets{ki}] = Upfun(CellSets{ki},CeSaSets{ki},TestSets{ki}(ii,:),Temp2);
            end            
            JG = [JG;Temp2-TestSets{ki}(ii,end)];
        end
    else
        [Cells,CeSas] = InitG(TrainData,mingen,t);
        DrawGrid(Cells,1)
        %plot(Tr1(:,1),Tr1(:,2),'b.')
        %plot(Tr2(:,1),Tr2(:,2),'k.')
        %plot(Tr3(:,1),Tr3(:,2),'g.')
        %plot(0.05,1.03,'b.')
        %plot(0.25,1.03,'k.')
        %plot(0.45,1.03,'g.')
        %text(0.06,1.03,'type 1')
        %text(0.26,1.03,'type 3')
        %text(0.46,1.03,'type 3')
        
        tt = size(TestData,1);
        for ii = 1 : tt
            %return affinity
            [Cells,CeSas,We] = CTestfun(Cells,CeSas,TestData(ii,:),dlttime);
            [Temp1,Temp2] = max(We);
            [Cells,CeSas] = Upfun(Cells,CeSas,TestData(ii,:),Temp2);
            JG = [JG;Temp2-TestData(ii,end)];
        end
        DrawGrid(Cells,2)
    end
    
    dgs = find(JG == 0);
    dr(j) = size(dgs,1) / tt * 100;
end
dr
Classification_accuracy = mean(dr)
Standard_deviation = std(dr,0,1) 


%%
%For high-dimension datasets
%Test stadge of continual learning
function [Grids Grs W] = CTestfun(Grids,Grs,data,dlttime)
%Data:n-dimension and m types
%data:[d1,d2,...,dn,time,label],time is the data acquisition
%Grids:[d1,d2,...,dn,p,q1,t1,q2,t2,...,qm,tm]
%di:dimension
%qi:number of samples in one cell
%ti:the last active time

di = size(data,2) - 2;
t = (size(Grids,2) - di - 1) / 2;
flag = 1;
ii = 1;

for ki = 1 : t
    Temp = find(Grids(:,di+2*ki+1)>0 & Grids(:,di+2*ki+1)<(data(end-1)-dlttime));
    if ~isempty(Temp)
        Grids(Temp(:,1),di+2*ki) = 0;
        Grids(Temp(:,1),di+2*ki+1) = 0;
    end
end

while flag
    [A] = IsSampleinGrid(Grids(ii,:),data);
    if ~isempty(A)
        mindj = max(Grids(:,di+1));
        if Grids(ii,di+1) < mindj  %High passage cell
            [TGrids,TGrs] = SuboneG(Grids(ii,:),data,mindj-1);
            ti = 1;
            flagT = 1;
            while flagT
                if ~isempty(TGrs{ti})
                    At = ti;
                    flagT = 0;
                else
                    ti = ti + 1;
                end
            end
            [adwgxh] = Fxlwgfun(TGrids(At,:),Grids,di);
            TE1 = Grids(adwgxh(:,1),di+2:2:di+2*t);
            dis = dist(Grids(adwgxh(:,1),1:di),data(1:di)');
            disi = TE1 ./ dis;
            if size(adwgxh,1) == 1
                W = disi;
            else
                W = sum(disi);
            end
            W = W ./ max(W);
        else
            [adwgxh] = Fxlwgfun(Grids(ii,:),Grids,di);
            XH = [adwgxh;ii];                              
            TE1 = Grids(XH(:,1),di+2:2:di+2*t);
            dis = dist(Grids(XH(:,1),1:di),data(1:di)');
            disi = TE1 ./ dis;
            W = sum(disi); 
            W = W ./ max(W);
        end
        flag = 0;
    end
    ii = ii + 1;
end
end


%%
%For high-dimension datasets
%Adjust the cell information
function [Grids,Grs] = Upfun(Grids,Grs,data,tp)
%Data:n-dimension and m types
%data:[d1,d2,...,dn,time,label],time is the data acquisition
%Grids:[d1,d2,...,dn,p,q1,t1,q2,t2,...,qm,tm]
%di:dimension
%qi:number of samples in one cell
%ti:the last active time

di = size(data,2) - 2;
flag = 1;
ii = 1;

while flag
    [A] = IsSampleinGrid(Grids(ii,:),data);
    if ~isempty(A)
        mindj = max(Grids(:,di+1));
        if Grids(ii,di+1) < mindj  %High passage cell
            [TGrids,TGrs] = SuboneG(Grids(ii,:),data,mindj-1);
            ti = 1;
            flagT = 1;
            while flagT
                if ~isempty(TGrs{ti})
                    At = ti;
                    flagT = 0;
                else
                    ti = ti + 1;
                end
            end
            [adwgxh] = Fxlwgfun(TGrids(At,:),Grids,di);
            smc = find(Grids(adwgxh(:,1),di+2*tp) > 0);
            if data(end-1) > min(Grids(adwgxh(smc(:,1),1),di+2*tp+1))
                Grids(adwgxh(smc(:,1),1),di+2*tp+1) = data(end-1);
                TGrids(At,di+2*tp) = 1;
                TGrids(At,di+2*tp+1) = data(end-1);
                Grids = [Grids;TGrids];
                Grs = [Grs;TGrs];
                Grids(ii,:) = [];
                Grs(ii,:) = [];
            end
        else
            [adwgxh] = Fxlwgfun(Grids(ii,:),Grids,di);
            smc = find(Grids(adwgxh(:,1),di+2*tp) > 0);
            XH = [adwgxh(smc(:,1),1);ii];
            if data(end-1) > min(Grids(XH(:,1),di+2*tp+1))
                Grids(XH(:,1),di+2*tp) = Grids(XH(:,1),di+2*tp) + 1;
                Grids(XH(:,1),di+2*tp+1) = data(end-1);
                Grs{ii} = [Grs{ii};data];
            end
        end
        flag = 0;
    end
    ii = ii + 1;
end
end


%%
%To find the adjacent memory cells 
function [xlwg] = Fxlwgfun(zx,Grids,di)
%Data:n-dimension and m types
%data:[d1,d2,...,dn,time,label],time is the data acquisition
%Grids:[d1,d2,...,dn,p,q1,t1,q2,t2,...,qm,tm]
%di:dimension
%qi:number of samples in one cell
%ti:the last active time

t = (size(Grids,2) - di - 1) / 2;
dqdj = zx(di+1);
dltyd = 0.5 / 2^dqdj;
TE1 = Grids(:,di+2:2:di+2*t)';
TE2 = find(TE1 < 0);
TE1(TE2(:)) = 0;
NTS = sum(TE1)';
AA = find(Grids(:,di+1) <= dqdj & NTS(:,1) > 0);
BB = Grids(AA(:,1),:);
tt = size(BB,1);
Gdzx = zx .* ones(tt,di+2*t+1);
Gdyd = 0.5 ./ 2.^BB(:,di+1);

fi = 0;
flag = 1;
while flag
    dltyd =(fi + 0.5) / 2^dqdj;
    zxyd = dltyd .* ones(tt,1);
    
    Gdab = abs(abs(Gdzx(:,1:di) - BB(:,1:di)) - Gdyd - zxyd);
    Gdmin = min(Gdab')';
    A = find(Gdmin <= 0.001 | Gdmin == 0);
    Grd = BB(A(:,1),:);
    Dis = dist(Grd(:,1:di),zx(1,1:di)');
    Diss = Dis - (di .* Gdyd(A(:,1),:) .* Gdyd(A(:,1),:)) .^ 0.5;
    dt = (di * dltyd * dltyd) ^ 0.5;
    B = find(abs(Diss - dltyd) <= 0.0001 | abs(Diss - dt) <= 0.0001 | Diss <= dt);
    
    if ~isempty(B)
        flag = 0;
    else
        fi = fi + 1;
    end 
end
%Return the sequence number of cell
xlwg = AA(A(B(:,1),:),:);
end

%%
%Divide the memory cell which has one sample to the smallest memory cell
function [TGrids,TGrs] = SuboneG(grid,data,xfdj)
%Data:n-dimension and m types
%data:[d1,d2,...,dn,time,label],time is the data acquisition
%Grids:[d1,d2,...,dn,p,q1,t1,q2,t2,...,qm,tm]
%di:dimension
%qi:number of samples in one cell
%ti:the last active time

di = size(data,2) - 2;
t = (size(grid,2) - di - 1) / 2;
dqdj = grid(di+1);
Temp = [];
Tempg = [];
TGrids = [];
TGrs = [];
if dqdj < xfdj
    dltyd = 1 / 2^(dqdj+1);
    yd = 0.5 * dltyd;
    for i = 1 : di
        k = 2^(i-1);
        Tei = [-yd.*ones(k,1); yd .* ones(k,1)];
        n = size(Tei,1);
        m = 2^di/n;
        Te = [];
        for j = 1 : m
            Te = [Te;Tei];
        end
        Temp(:,i) = grid(i) + Te;
    end    
    Temp = [Temp (dqdj+1).*ones(2^di,1) zeros(2^di,2*t)];  
    for i = 1 : 2^di
        Tempg{i,1} = [];
    end
    flagt = 1;
    while flagt && j <= 2^di
        [A] = IsSampleinGrid(Temp(j,:),data);
        if ~isempty(A)
            Tempg{j} = data;
            Temp(j,di+(2*data(end))) = 1;
            Temp(j,di+1+(2*data(end))) = data(end-1);
            [TGrids,TGrs] = SuboneG(Temp(j,:),data,xfdj);
            Temp(j,di+1) = 999999;
            flagt = 0;
        else
            j = j + 1;
        end
    end
    TGrids = [TGrids;Temp];
    TGrs = [TGrs;Tempg];
else
    dltyd = 1 / 2^(dqdj+1);
    yd = 0.5 * dltyd;
    for i = 1 : di
        k = 2^(i-1);
        Tei = [-yd.*ones(k,1); yd.*ones(k,1)];
        n = size(Tei,1);
        m = 2^di/n;
        Te = [];
        for j = 1 : m
            Te = [Te;Tei];
        end
        Temp(:,i) = grid(i) + Te;
    end    
    Temp = [Temp (dqdj+1).*ones(2^di,1) zeros(2^di,2*t)];   
    for i = 1 : 2^di
        Tempg{i,1} = [];
    end
    j = 1;
    flagt = 1;
    while flagt && j <= 2^di
        [A] = IsSampleinGrid(Temp(j,:),data);
        if ~isempty(A)
            Tempg{j} = data;
            Temp(j,di+(2*data(end))) = 1;
            Temp(j,di+1+(2*data(end))) = data(end-1);
            flagt = 0;
        else
            j = j + 1;
        end
    end
    TGrids = [TGrids;Temp];
    TGrs = [TGrs;Tempg];
end
AA = find(TGrids(:,di+1) == 999999);
if ~isempty(AA)
    TGrids(AA(:,1),:) = [];
    TGrs(AA(:,1),:) = [];
end
end


%%
%initializing memory cells
function [Grids,Grs] = InitG(Data,mingen,t)
%Grs records the samples in every memory cell, it is used to check the result of Grids and continual learning
%Data:n-dimension and m types
%data:[d1,d2,...,dn,time,label]
%Grids:[d1,d2,...,dn,p,q1,t1,q2,t2,...,qm,tm]
%di:dimension
%p:passage
%q:number of samples in one cell;qi type i
%t:the last active time
di = size(Data,2) - 2;                                                     %The dimensionality of Data

GC = 0.5 .* ones(1,di);
%1-di,center of grid; di+1,passage, di+2,sample size of type1, di+3,sample size of type2
Grids(1,:) = [GC zeros(1,2*t+1)];                                            %The initial center, passage, number of type1, number of type2
Grs{1,1} = [];

%Adaptive to generate memory cells
[Grids,Grs] = AdjustG(Grids,Grs,0,Data,mingen);                            %Adaptive to generate memory cells
%DrawGrid(Grids,1)
end


%%
%Adjusting memory cells and samples
function [Grids,Grs] = AdjustG(Grids,Grs,gdj,Data,mingen)
%gdj is the passage of memory cell for adjusting
%Data:n-dimension and m types
%data:[d1,d2,...,dn,label]
%Grids:[d1,d2,...,dn,p,t1,t2,...,tm]
%di:dimension
%t:types
tt = size(Grids,1);
di = size(Data,2) - 2;
for ii = 1:tt
    [A] = IsSampleinGrid(Grids(ii,:),Data);
    Grs{ii} = A;   
    if ~isempty(A)
        pas = Grids(ii,di+1);
        N = unique(A(:,end));
        if pas < mingen            
            [TempG,TempS] = Xfwgfun(Grids(ii,:),gdj,Grs{ii},mingen);              %Divede memory cell
            Grids(ii,di+1) = 999999;
            Grids = [Grids;TempG];
            Grs = [Grs;TempS];
        else
            tj = size(N,1);
            for jj = 1 : tj
                TA = find(A(:,end) == N(jj));
                Grids(ii,di+2*N(jj)) = size(TA,1);
                Grids(ii,di+2*N(jj)+1) = A(TA(1),end-1);
            end
        end
    end
end
AA = find(Grids(:,di+1) == 999999);
if ~isempty(AA)
    Grids(AA(:,1),:) = [];
    Grs(AA(:,1),:) = [];
end
end


%%
%The sample in or not in a memory cell
function [data] = IsSampleinGrid(grid,Data)
%Data:n-dimension and m types
%data:[d1,d2,...,dn,label]
%Grids:[d1,d2,...,dn,p,t1,t2,...,tm]
%di:dimension
%t:types
di = size(Data,2) - 2;
dltyd = 1 / 2^grid(di+1);
flag = 1;
i = 1;
while flag && i <= di
    L = grid(i) - 0.5 * dltyd;
    R = grid(i) + 0.5 * dltyd;
    A = find(Data(:,i) <= R & Data(:,i) >= L);
    if isempty(A)
        flag = 0;
    else
        Data = Data(A(:,1),:);
    end
    i = i + 1;
end

if flag
    data = Data;
else
    data = [];
end

end


%%
%Divide memory cell
function [Grids,Grs]= Xfwgfun(grid,xfdj,data,mingen)
%Data:n-dimension and m types
%data:[d1,d2,...,dn,label]
%Grids:[d1,d2,...,dn,p,t1,t2,...,tm]
%di:dimension
%t:types
dltyd = 1 / 2^(xfdj+1);                                                    %The length of other divide memory cell
yd = 0.5 * dltyd;
di = size(data,2) - 2;
t = (size(grid,2) - di - 1) / 2;
for i = 1 : di
    k = 2^(i-1);
    Tei = [-yd.*ones(k,1); yd.*ones(k,1)];
    n = size(Tei,1);
    m = 2^di/n;
    Te = [];
    for j = 1 : m
        Te = [Te;Tei];
    end
    Temp(:,i) = grid(i) + Te;
end
Temp = [Temp (xfdj+1).*ones(2^di,1) zeros(2^di,2*t)];
for i = 1 : 2^di
    Grs{i,1} = [];
end
[Grids,Grs] = AdjustG(Temp,Grs,xfdj+1,data,mingen);
end


%%
%Drawing memory cell in [0,1]^2
function DrawGrid(Grids,k)
%k is the figure number
[tt1 tt2] = size(Grids);
tt = tt2 - 3;
kt = ceil(tt / 6);
cl = 1/(kt+1);
mt = max(Grids(:,end))+1;

figure(k)
hold on
for fi = 1 : tt1
    dltyd = 1 / 2^Grids(fi,3);
    L = Grids(fi,1) - 0.5 * dltyd;
    R = Grids(fi,1) + 0.5 * dltyd;
    U = Grids(fi,2) + 0.5 * dltyd;
    B = Grids(fi,2) - 0.5 * dltyd;
    for kk = 1 : 2 : tt
        TE1 = Grids(fi,3+1:2:3+tt);
        TE2 = find(TE1>0);
        TET = sum(TE1(TE2(:)));
        if Grids(fi,3+kk) > 0 && Grids(fi,3+kk) == TET %one type
            cc = [0 0 0];
            ft = floor(kk/2) + 1;
            cc(mod(ft-1,3)+1) = cl * (floor((ft-1)/3)+1) * (Grids(fi,3+kk+1)+1)/mt;
            X = [L L R R];
            Y = [B U U B];
            fill(X,Y,cc)
        end
        if Grids(fi,3+kk) > 0 && Grids(fi,3+kk) < TET %more than one type
            cc = [0.25 0.5 0.75];
            X = [L L R R];
            Y = [B U U B];
            fill(X,Y,cc)
        end        
    end
    plot([L;L;R;R;L],[B;U;U;B;B],'g-');
end
plot([0,1],[1,1],'k-');
plot([0,1],[0,0],'k-');
plot([0,0],[1,0],'k-');
plot([1,1],[1,0],'k-');
axis equal;
axis off
set(gcf,'color','w')
end


%%
%normalize data to [0,1]n
function [Dn] = norma(D)
if max(D) == 0
    Dn = zeros(size(D,1),1);
else
    Dn = (D - min(D)) ./ (max(D) - min(D));
end
end


%%
%normalize data to [0.0001,0.9999]n
function [Dn] = normat(D)
if max(D) == 0
    Dn = zeros(size(D,1),1);
else
    Dn = 0.0001+0.9998.*(D - min(D)) ./ (max(D) - min(D));
end
end